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Abstract 



We investigate nuclear matter on a cubic lattice. An exact thermal formal- 



Q^ ■ ism is applied to nucleons with a Hamiltonian that accommodates on-site and 



next-neighbor parts of the central, spin- and isospin-exchange interactions. 
We describe the nuclear matter Monte Carlo methods which contain elements 
from shell model Monte Carlo methods and from numerical simulations of 



S^ . the Hubbard model. We show that energy and basic saturation properties of 

nuclear matter can be reproduced. Evidence of a first-order phase transition 
from an uncorrelated Fermi gas to a clustered system is observed by comput- 
ing mechanical and thermodynamical quantities such as compressibility, heat 
capacity, entropy and grand potential. We compare symmetry energy and 
first sound velocities with literature and find reasonable agreement. 

I. INTRODUCTION 

Properties of nuclear matter have been deduced from different approaches. The volume 
term of the semi-empirical mass formula [|l],H predicts a binding energy of 16 MeV, and 



calculations of finite nuclei estimate the equilibrium density to be po = 0.16 fm~^. However, 
properties of finite nuclei are strongly influenced by finite size effects like the surface effect, 
and it is therefore difficult to estimate energies and saturation density of nuclear matter from 
nuclei. Many-body calculations of nuclear matter are based on sophisticated potentials (such 
as the Argonne AV14 and AV18 and Urbana UV14 potentials) and use Bethe-Brueckner- 
Goldstone theory l^-^f and hypernetted chain approximations [^0 to calculate ground state 
properties. Lattice gas calculations |]8 HlO|] attempt a thermal description of nuclear mat- 
ter. They work with much simpler Hamiltonians, incorporating isospin-1 or Hubbard-like 
interactions. These calculations aim at the investigation of a liquid-gas phase transition of 
nuclear matter expected to take place at subnuclear densities and low temperatures. They 
are classical, not quantum mechanical, putting in kinetic terms by hand or sampling them 
from a Maxwell-Boltzmann distribution. These types of calculations use Monte-Carlo-like 
algorithms and show that the inclusion of a kinetic term is crucial to observe a phase tran- 
sition. 

This paper describes first results of a calculation of infinite nuclear matter that com- 
bines both the usage of a more realistic Hamiltonian and the exact, thermal treatment of 
the many-body problem on a lattice. In the last few years, the shell model Monte Carlo 
(SMMC) method has been successfully developed |]TT1-|T5[] to give a powerful alternative 



to direct diagonalization procedures which suffer from the fact that the many-body space 
scales so unfavorably with the number of single-body states considered. Direct diagonaliza- 
tion methods can only address very light nuclei or nuclei with a closed shell and only a few 
valence nucleons. The SMMC avoids this combinatorial scaling (in storage and computation 
time) and makes it possible to investigate structural properties of nuclei far beyond the few- 
nucleon system. The SMMC enforces the Pauli-principle exactly, and concentrates on the 
evaluation of thermal averages of observables. This would be the main purpose of a nuclear 
matter investigation too: Not focusing on obtaining a wave function, a thermal formalism is 
useful for a study of nuclear matter, because the equation of state is of main interest, which 
clearly depends on density and temperature. Moreover, the consideration of a large piece 



of infinite nuclear matter in coordinate space reduces finite size effects that appear after 
imposing periodic boundary conditions. A formalism written in momentum space has the 
disadvantage that two- or many-body correlations cannot be calculated directly: Clustering 
(and therefore a possible liquid-gas transition) is not as easily calculated and observed as in 
the coordinate space representation. 

The following concept is pursued for our nuclear matter calculation: The quantum me- 
chanical and exact treatment with the full Hamiltonian, kinetic and potential term, should 
be a prerequisite for a successful description of the physical system. In a coordinate repre- 
sentation nucleons shall interact with a potential that eventually comes as close to a realistic 
nuclear interaction (like AV18) as possible. The partition function along with observables of 
interest shall be calculated in the grand canonical ensemble, in order to control temperature 
T and density p. The latter is to be adjusted on average via the chemical potential /i. The 
many-body problem shall be solved exactly using Monte Carlo methods similar to those 
used in the SMMC applications. At the same time, realizing that the emerging equations 
eventually have to be solved on a computer, one should take into account that space will 
be discretized, and advantage should be taken of the available technology that has been 
employed for the Hubbard and other models in condensed matter physics. This paper is to 
be viewed as a first step of a full thermal description of nuclear matter in which we con- 
strain our potential parameters to a reasonable shape of the energy as a function of density, 
including the correct saturation point. 

II. THEORY OF NUCLEONIC MATTER ON A LATTICE 

The general concept of the nuclear matter calculation consists of nucleons interacting via 
a variety of components of the nuclear two-body potential. While it should be the ultimate 



goal to use a potential that fits the nucleon-nucleon scattering data best |I6[, at the first 
stage we concentrate on few parts of the interactions, namely central, spin- and isospin- 
exchange. The degrees of freedom of the nucleon are its spin, isospin as well as the spatial 



coordinate. 

Subnuclear degrees of freedom are not explicitly incorporated. The lightest meson, the 
pion, facilitates an interaction with a range of r ^ 1.4 fm which is of same order as the lattice 
spacing of the applications in this paper. Since the system is ultimately regularized on a 
lattice, the argument can be made that all subnuclear degrees of freedom are integrated out, 
resulting in a strong on-site and weaker next-neighbor interaction. The lattice spacing, here 
an additional fitting parameter like the potential parameters, is chosen to be of a = 1.842 fm. 
This particular lattice spacing sets the half- filling of the lattice at p = 2po = 0.32 fm"^. 
Other settings of fillings have been tried, but turned out to reproduce the saturation curve 
less well. 

In this section we specify the Hamiltonian of the system and describe the nuclear matter 
Monte Carlo method (called NMMC hereafter), which consists of the thermal formalism 
to express the grand canonical partition function as an integral over single-body evolution 
operators. At its center stands the Hubbard-Stratonovitch transformation, which is used 
to reduce the many-body problem to an effective one-body problem. The details of the 
Monte Carlo procedure, which is used to evaluate the resulting multi-dimensional integral, 
are explained. 

A. Hamiltonian 

We consider a three-dimensional cubic lattice of spacing a and assume periodic boundary 
conditions, which result in a three-dimensional toroidal configuration. The coordinate x and 
the momentum p are discretized as 

X — > am = Xm, (1) 



'^^(l^)*"P" P' 



such that 



- - 27r . 

X ■ p = — X integer, (3) 

where N is the number of lattice points in each spatial direction, and fh and k are vectors 
with integer components. 

The nucleons have mass miy, spin o" = ±| and isospin r = ±| . The Hamiltonian, 

n = t + v, (4) 

is expressed in second quantization and contains kinetic and potential operators. The kinetic 
term is written as 

/C = --^$: /df C(^)VV..(:?). (5) 

The fermion operator ■?/'^^(x) creates a nucleon of spin and isospin (a, r) at location x, while 
its adjoint iparix) destroys it. This equation is discretized on the lattice by the symmetric 
3-point formula for the second derivative, and the integral is replaced by a finite sum, which 
results in 

iC = -tQ^a^ ^ ^l^{Xn)[i'<TT{Xn + a'ei)-2^„r{Xn)+i'aT{Xn- aci)] (6) 

err ^„,i=l---3 



with 



Here, the orthogonal unit vectors {cj} span the three-dimensional space. 

While the form of the nuclear potential is generally given, we here are limited by current 
computational constraints. The treatment of a full Hamiltonian, as it is represented in 
the Argonne potential, for example, is computationally impossible with currently available 
computer power, but may be feasible in a few years. We chose 

V = Vc + %. (8) 

The first part is the central potential (Vc); followed by the spin-exchange (Vo-)- The general 
form for the scalar potential. 



can be written in terms of the density 

P(^) = E P'^rix) = E i^trix)^<^r{x). (lO) 



ar cfT 



The purpose of doing so is to cast the potential in hnear and quadratic terms, as the 
Hubbard-Stratonovitch transformation can only be performed on quadratic terms. Using 
the fermion anticommutation relation, the potential then becomes 

Vc = ^ / df I df V,ix - ^)p{x)p{x') -\j dx V,{0)p{x). (11) 

The last term is the self-energy and is a consequence of the Pauli principle. The discretized 
version of this equation is 

^c = Y E VA - <)/5(f„)p(f;) - y E ^c(0)p(f„). (12) 

We assume a Skyrme-like on-site and next-neighbor interaction 

K (f„ - x'^) = V}'^5 (f„ - x'J + V;(2) {vy (f„ - <)) , (13) 

whose discretized form is 

Vc{Xn - X'J = -^6s„,£'^ + -^ E {^x„+ae,,x'„ " 25^„,,f^ + (5^„_ae„x'4 • (14) 



The parentheses in Eq. (|T3|) indicate that the Laplace operator only acts on the (5-function, 
but not on any following parts. Inserting equation (|T^) into (0) gives 

^- = ^ E a^KXnf - ^Y~ E E (P (^" + «e;) - P {Xn)) 

Xn Xn 1 J- 

-l(yc'''-Q^)T.Pi^n), (15) 

Here, we applied periodic boundary conditions. 

The spin-exchange part of the potential is handled in a very similar way. Starting from 

^<^ = 9 E / df / df V^J^(x)V^J,^,(f' )Kr(^- ^)%kA ■ ff^'T'K'y'ipt,'X'{x')ipKx{x), (16) 



kAk'A' 



we write the potential in form of spin densities 



PT\^) = E %(^)4riA^-A(x), a = 0, +, -, 



(17) 



^tkA 



.(") 



where (tL^;^ are the elements of a generalized Pauli spin-isospin matrix, which acts on a 
4-vector representing all spin/isospin states of a nucleon. We assume the the same spatial 
dependence ([T3|) as for the central part, and finally arrive at 

^- = ^ E «' (Pi°^ i^nf + 2 (p(+) (f„) + /5(-) (f„))') 

- ^ E E ({p'^H^n + ael) - pW(x„))' 

+ 2 (p^+) (f„ + aci) + pi-^ (xn + aci) - p^+) (f„) - p^^^ (fn))^) 

-^(K^°^-6^)EP(^^n)- (18) 

Other components of the potential can be included in a similar way. 



B. Nuclear Matter Monte Carlo Method 



In order to study thermal properties of nuclear matter, the grand canonical partition 
function at a given temperature T = j3~^ needs to be determined. 



Tr 



e^Y>[-l3\n-Y.PrKA 



Tr 



u 



(19) 



with A/o-T = Sxn i^tri^nj'ipcTTixn) and fir &§ the isospin-dependent chemical potential. U is 
called the imaginary-time evolution operator of the system and is a many-body operator. 
In the present study the Hamiltonian 7i contains one- and two-body operators as described 



Section [II A| , and the trace is taken over all many-body states as indicated by a caret. The 
partition function Z is an exponential over all one- and two-body operators (and therefore 
interactions) present in the system. It is impossible to deal with the partition function Z 
in this form, because the number of many-body correlations that have to be kept track of 
grows rapidly with system size. We therefore seek an expression for Z that is based on 



a single-particle representation, and we will replace the many-body problem with that of 
non-interacting nucleons that are coupled to a heat bath of auxiliary fields. This involves 
rewriting Z as a multi-dimensional integral. 

We start by dividing the evolution operator into rtt time slices: 



U = exp (-P (n~J2 /^-^-) ) = exp (-Af3 (n-J2 f^rAr) 



(20) 



with A(3nt = j3. The Trotter approximation |T7|JT8|] is used to separate one-body (kinetic 
energy and chemical potential) in 5 = Ai — '2Za,T y^T^m and two-body terms (potential) in 

exp I -A/5 [n-Y. [iMar\ j = exp (-A/3 (s ^V)) ^ exp (-A/?5) exp (-A/3V) . (21) 



Equation ( |2lD is valid to order A/5, but becomes exact in the limit A/5 -^ 0. 

The propagator of each time slice for the potential, exp(— A/3V), is manipulated by 
applying the Hubbard-Stratonovitch (HS) transformation p!9| , |20| on each term, replacing it 
with a multi-dimensional integral over a set of auxiliary fields. 

As an example, we describe the transformation by taking the on-site part of Vc at one 
particular site x^. Using a = A/3V"J°V2 and defining 



5n 



±i if a > 
±1 if a <0, 
the propagator for this single interaction is written as 



(22) 



MJiXr, 



exp(-A/5^p2(^^ 




TT 



2 
dxexp(- 



— / dxexp(— ap [x^ 

TX J-oo 



\a\ 



X^ + '2SaXP {Xra) 



\a\{x + SaP{Xm)f 

(23) 



The last line of Eq. ( p3| ) reveals that the evolution operator is now expressed in terms of 
the exponential of a one-body operator and an integration over the auxiliary field x- It 
has become a one-body propagator that corresponds to non-interacting nucleons coupled to 
this field. Since the integral is calculated with Monte Carlo methods, the field fluctuates 
according to a weight that is to be specified, hence the picture of a heat bath. 
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It has to be emphasized that p (x^) here represents a one-body operator for a subset (Vc 



in this case) of the full interaction. Each quadratic term in ([15|) and (|T8D has to be replaced 
by an integral. At a given lattice site Xm-, there are twelve auxiliary fields to form the full 
interaction, four for V^ and eight for V„. Nucleons are now coupled to a big ensemble of 
auxiliary fields through which the interaction of the nucleons is mediated. The Af/'s are then 
multiplied together to form the evolution operator for one time slice U {f3m) at (3m = mAf], 
which is expressed only in terms of single-body matrices, and ultimately, all time slices are 
multiplied together to form U: 



U 



exp {-Af3n)Y' = J^ix]G ix) U^ (/?, 0) (24) 



with the integration measure 



'la 



^[x]=UUU'iXm,,.A^- (25) 



m=l x„ i ' 



The Ui = ApVi/2, Vi G Vlp\ V^i^\ '^!P\ ' ' ' ? are the interaction-specific coupling strengths 
of auxiliary fields to nucleons, and the index i enumerates all fields at a particular site. The 
Gaussian factor G is given by 

nt 

G{X) = n nnexp(-|a.|xl.^„,) , (26) 

m=l Xn i 

and the one-body propagator is 

U^ iP, 0) = U i(3„,) If (/?„,_i) ■■■U (A) . (27) 

Note that Eq. (p^) only becomes exact in the limit of an infinite number of time slices, 
rif —>■ oo. For a finite rit, the Hubbard-Stratonovitch approximation is valid only to order 
A/3. 

In the practical implementation of Eq. (p3D, a discrete Hubbard-Stratonovitch trans- 
formation is used instead of the continuous form because it turns out to use much fewer 
de-correlation sweeps, as explained below. 



A thermal observable (O) is expressed as [ll2| , ^ 



(O) 



:Tr 



dexp(-l3(H-J2fIrJ^^r]] 



Iv[x]Gix){dixmix) 
J-DlxMxKix) 



(28) 



and has the integration measure of Eq. (p5|) and Gaussian factor of Eq. p6|). The expecta- 
tion value of any operator in second quantization can be obtained with the help of Wick's 



theorem, and the resulting one-body densities are |12 



(C (^n) ^.v (x™))x = { [1 + U^ iP, 0)]-' U^ if3, 0)} 



(29) 



{cr'T',Xm),io-T,Xn) 

The bold face U^ (/^, 0) is the single-body matrix representation of U^ (/5,0). Observables 
of the system are chosen to be the number of neutrons and protons and all components of 
the energy. 

The integrals in ( pS] ) are evaluated using the Metropolis algorithm |2^. The basic idea 
involves sampling the integrand of 



{0{X)) 



Tr 



6f/^(/?,o) 



Tr 



U^ i/3, 0) 



(30) 



within the boundaries of the integration volume according to a positive-definite weight 



Wix) 



\G{x)^{x)\ for continuous HS 



\ax)\ 



for discrete HS, 



(31) 



with 



e(x) = Tr f/J/5,0) =det[l + U^(/3,0)]. 



(32) 

The last equality can be proven by expanding the determinant |^ . Samples are taken by a 
random walker that travels through x-space, taking a new value Xnew from the previous one 
Xoid if the ratio 

WiXnew) 



W{xom) 



(33) 



is larger than one, or else, if r < 1, taking on Xnow with probability r. It can be shown ||2^ 
that the sequence of values the walker takes is distributed according to the weight function 
W{x), which is typically chosen to be as close to the integrand as possible to increase the 
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efficiency of tlie procedure. Since tlie consecutive steps are correlated, tlie walker has to 
travel several steps before another sample is taken to de-correlate them. The average of an 
observable (^) is then simply 

id) - ^^ (34) 

in terms of its ith Monte Carlo sample {0)i and 

I ^Wf for continuous HS 
$= J>^^ (35) 

r^^ for discrete HS. 

Note that $, which is just the sign of the weight W, can generally be negative or even 
complex. 

The numerical determination of Eq. (0), which is the Monte Carlo equivalent of the in- 
tegrals in Eq. (|28D , can be difficult in certain situations, even with Monte Carlo methods: If 
Sa = ±^ (which generally corresponds to a repulsive on-site and an attractive next-neighbor 
interaction, cf. Eq. (p2D ), propagators for the potential (Eq. (p3D ) contribute negative or 
complex elements to U^{j3,0) (see Eq. (|3^)). The integrands in both numerator and de- 
nominator are oscillatory, and the integrals can add up to small numbers. A numerical 
evaluation with Monte Carlo methods causes large uncertainties because the methods are of 
a stochastic nature, and the number of samples in a computation remains finite. This is a 
complication associated with these methods when the Hubbard-Stratonovitch transforma- 
tion is used. It is referred to as the Monte Carlo sign problem. A pragmatic solution has 
been used for the shell model Monte Carlo method to handle this complication [|15 . 



There has been significant effort in stabilizing and optimizing the Metropolis algorithm 
for lattice calculations, as they have been heavily used for models of interacting electrons in 
condensed matter physics. Many of the techniques have directly been applied to NMMC, 
because the models are similar. Besides using the checkerboard breakup |21]] technique for 
kinetic and spin-exchange parts, we use the Green's function algorithm described in [^ to 
reduce the computational burden. 
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III. NUMERICAL RESULTS 

We now show that for symmetric nuclear matter (SNM) the energy per particle can be 
reproduced quite well over a wide range of densities, and the energy for pure neutron matter 
(PNM) for the same potential is discussed. Then, several observations are presented that 
give evidence of a first-order phase transition from a Fermi gas to a clustered system at a 
critical temperature Tc ~ 15 MeV. Furthermore, the symmetry energy and first sound are 
discussed. 

The extensive search in the space of potential parameters included all components of the 
central part and spin-exchange. The effort focused on reproducing saturation density and 
energy correctly. The project is considered to be a first step towards a realistic calculation 
as indicated earlier. A more realistic calculation has to contain more parts of the nuclear 
potential and the spatial resolution has to improve; a perfect fit over a wide range of densities 
should not be expected at this point. The fit has been performed in such a manner that 
the saturation energy and density is reproduced, and that the overall energy curve has a 
reasonable form: for sub-saturation densities, the matter should be unstable, while for p > po 
SNM should evolve in an unbound state {E/A > MeV for p > 0.4 fm~^). 

The sign problem unfortunately forces the use of a nuclear potential that might contradict 
the usual physical understanding and intuition based on few-nucleon potential models. It 
is generally known that the central potential has a strong repulsion for short distances and 
features a long-range attraction. Here, the desire to avoid the sign problem generates the 
opposite: on-site attraction and next-neighbor repulsion. On the other hand, an on-site 
attraction and next-neighbor repulsion may not be unreasonable given the fact that there 
have been several mean-field calculations of nuclear matter with the Skyrme forces. Skyrme 
forces simulate the interaction with a 6-like attraction and a V^5-like repulsion. In the 
lattice discretization of this investigation, the position of the nucleons at the same site is 
only determined up to a cube of size a. Therefore, the on-site potential parameter can be 
seen as an average potential within that cube, and by tuning the lattice spacing accordingly, 
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it could be possible that this parameter is negative. The exact definition of the parameter 
depends on a regularization scheme. In such a scheme the Schrodinger equation has to 
be solved on a lattice, and by identifying scattering amplitudes, one could determine the 
potential parameters from scattering lengths and effective range |Q. With respect to the 
lattice spacing a two approaches can be taken: First, we constrain ourselves to a description 
with a fixed number of lattice sites. Then a becomes a free fitting parameter like the potential 
parameters, and the latter would have to be interpreted as an average potential, as noted 
above. In the second approach, a is an discretization parameter for the potential, and the 
ultimate goal would be to increase the number of lattice points with decreasing a, getting 
a smooth parameterization of the potential. If, in that case, a positive on-site parameter 
is used, emulating a hard core repulsion, one has to deal with oscillatory integrands and 
commensurately large error bars. 

The following potential parameters were obtained: 

1/(0) = -181.5 MeVfm^ (36) 

1/(2) = 37.8 MeVfm^ (37) 

1/(0) = -31.25 MeVfm^, (38) 

V;(2) = O.OMeVfml (39) 

All calculations were done with this set of parameters. The lattice has a spacing of 

a = 1.842 fm, (40) 

tuned such that quarter filling of the lattice is at saturation density po = 0.16 fm~^. In 
this paper, the lattice spacing is an additional fitting parameter, and several other settings 
have been tested. However, quarter- and half-filling at saturation density have a special 
significance because certain lattice occupations of the nucleons result in an energetically 
favored configuration. 

Because of limited CPU time, the calculation is restricted to 4 x 4 x 4 lattices for the 
moment. This lattice comprises 10^^ many-body states, and 11520 auxiliary fields are used 
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for this set of parameters. All calculations are prepared by a pre-thernialization of the 
system for 100 steps before we took measurement samples. Between measurement samples, 
15 de-correlation steps have been taken to guarantee statistical independence of the samples. 



The autocorrelation of k consecutive samples [23 



Co {k) - -^aW^rja^^ (4i) 

with i being the summation index over samples, has been monitored for all observables O 
and was held below 10%. 

We are also restricted by the fact that the Monte Carlo simulations cannot be extended to 
arbitrarily low temperatures. Even though the numerical routines are quite stable, it is not 
possible to add an arbitrary number of time slices, since it involves more and more matrix 
multiplications which become increasingly numerically unstable. In the present case we take 
AP = 0.01 MeV-^ and were able to go down to a value of /? = 30 x A/? = 0.3 MeV"^ without 
running into numerical instabilities. Thus, the temperature range of the investigation is 

3.0 MeV < T < 100 MeV. (42) 

Fig. |I] shows the best fit we obtained. With decreasing temperature, the system develops 
a minimum at p = 0.32 fm~'^ first, which is most pronounced between 10 — 14 MeV, before 
it shifts to lower densities. At T = 3.3 MeV and T = 5.9 MeV the minimum is very 
broad, making matter softer (see also compressibility. Fig. ^ below). For high temperatures 
and/or high density, the simulation suffers from the fact that it runs out of model space: 
At T = 50 MeV the system behaves almost like a Fermi gas and the energy per particle 
should behave like ~ p^/^. Yet, the curve bends down. Also, for all other temperatures, 
the curves converges to the energy of the full lattice state, E/A = 5.96 MeV, as density 
increases. For sub-saturation densities the model gives more binding if compared to other 
calculations (see, for example, Refs. and [|^), and the energy is not as high for densities 
beyond saturation. At p = 0.32 fm~^, E/A as a function of temperature has a minimum 
at T ^ 10 MeV which means that at even lower temperatures E/A increases again. This 
contradicts intuition because it would mean that the system is in an unphysical state. 
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The last issue needs further explanation. The energy per particle is not the correct 
quantity in order to address the question of stability. Particles fluctuate in and out of 
the system differently at different temperature, leaving the average number of particles 
unchanged, but contributing to the two-body part of the Hamiltonian. If, however, the 
grand potential is plotted (see Fig. 0), 

f](/3,/i) = -rinZ(/?,/x), (43) 

with 

In Z (/?, /x) - In Z (0, /i) = - /'' d(3'E (/?', /x) , (44) 



it turns out to actually be a monotonic function of temperature, with a slight deviation at 
/i = 11.0 MeV where the negative slope of f2, the entropy 

becomes zero between 10 MeV and 14 MeV and positive again for even lower temperatures. 
This is a slight anomaly (see also Fig. |^) which may have been caused by the onset of 
numerical instabilities at low temperatures or the fact that the lattice spacing is so big and 
the number of sites so small that the discretization of space is not accurate enough. 

We now introduce several observations which indicate that the system may undergo a 
first-order phase transition towards a clustered system when the temperature is lowered. 
First, we investigate changes in density with respect to the chemical potential /i. It is well 
known that they are proportional to particle fluctuations 



2 ^ d{N) 






(46) 

T,V 



T,V 

Such fluctuations are typical for first-order phase transitions and indicate that particles move 
between the two phases without energy cost. For an infinite system, the fluctuations should 
diverge, but not for a finite system. In the present case, we expect particles building clusters 
and breaking them up again, so one phase — the gas phase — would be that of independent 
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particles, the other one that of clusters. Since we observe the single particle density, a% 
describes the fluctuations in the gas phase in which p is linear in p. At T = 100 MeV, we 
are in the gas phase with nucleons behaving like a Fermi gas. Therefore, to simplify the 
graphs, we have flrst fltted the data at T = 100 MeV to a linear function, 

Pfit = Oflt + &flt X /i, (47) 

and then subtracted this function from all data points of all temperatures, deflning a function 
of temperature and chemical potential 

/(T,/i) = p(r,/i)-pfit. (48) 

The upper panel of Fig. | shows the outcome of this procedure. We then take the derivative 
of / (T, /i) with respect to /i and multiply with T, and this is shown in the lower panel of 
Fig. 0. The fluctuations show a pronounced maximum for T = 14.3 MeV and /i ~ —8 MeV, 
while they are low for T = 3.3 MeV and T = 100 MeV. The phase transition seems to occur 
somewhere between T = 8 MeV and T = 20 MeV. 

Another quantity that suggests the existence of a transition is the compressibility which 
is given by 

2 d^E/A 



K = 9p' 



(49) 

P=Psat 



dp"^ 

where psat is the saturation density. We have fltted the minima of each energy curve to a 
quadratic function 



E 



(^p) =a + bx {p- psat)^ , (50) 

fit 



and determined the compressibility as k = ^8p'^^^ x b. All data points in Fig. ^were obtained 
with a x^ per degree of freedom of less than 1.5. Again, a maximum in compressibility (which 
is in fact an incompressibility) is observed at T ^ 14 MeV: The clusters that form repel each 
other through the next-neighbor interaction which is repulsive. At T < 14 MeV, matter 
becomes softer again due to a broadening of the minima in E/A. This can be explained if 
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one assumes that the system becomes more dilute. Note that the values of psat change with 
temperature. 

Finally, we present the heat capacity and entropy of the system. For a first-order phase 
transition, the continuous and infinite system shows a divergence in the heat capacity 

dE 



'^- ar 



(51) 

V 



and a discontinuity for the entropy (with an infinite derivative at Tc), 

S{P,^i)=\nZ + [3{H{[3,^^))-[3^^{N{P,^^)). (52) 

For a finite system only a maximum in the heat capacity is expected, and a relatively sharp 
drop in entropy with decreasing temperature. Both facts can be verified in Fig. ^ The 
heat capacity suggests a critical temperature of Tc = 15 MeV, as does the entropy. For the 
graphs of entropy one has to keep in mind that the system investigated is quite small (it 
is a 4 X 4 X 4 lattice only), but two levels, 5 ^ 75 from T = 5 MeV to T = 12 MeV and 
S* = 175 from T = 20 MeV to T = 30 MeV with a steep decrease in between, can definitely 
be identified. To show that this is indeed a first-order phase transition, the calculation has 
to be repeated for a larger number of lattice sites, and then it has to be demonstrated that 
the drop between the two levels becomes steeper and steeper, finally resulting in a step-like 
function. This will have to be left for a future project. Studies of a lattice gas model ||2^ have 
shown that finite size effects do induce anomalies in physical quantities, and a first-order 
phase transition rather appears as a second-order one. The anomaly below T = 10 MeV has 
been addressed when discussing the grand potential. However, the latter quantity shows a 
qualitative behaviour at /i = 11 MeV that is expected for a phase transition (cf Fig. |^). The 
infinite system has a kink (the derivative is not continuous) in the grand potential at the 
critical temperature. Consequently, all quantities consistently suggest a phase transition at 
a critical temperature of Tc ~ 15 MeV. 

In Fig. ^ we show the energy per particle for pure neutron matter. The uncertainties 
for this case are much larger than for symmetrical nuclear matter. As a potential, we 
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used the parameters obtained from the fit to symmetric nuclear matter, even though we 
could have fitted the potential parameters for this case anew, including an isospin-exchange 
potential. Therefore we view the results for pure neutron matter more as a test to see how 
well the given potential already reproduces the energy. Note that the slopes of the curves at 
high temperatures are not negative as it is for symmetrical nuclear matter. But clearly, we 
cannot conclude that the energies at T = 3.3 MeV have converged to that of the ground state 
because the curve differs quite a bit from that of T = 5.9 MeV. At the lowest temperature 
they are 4 — 5 MeV higher than those of the ground state as calculated in Ref . |^ , but the 
general shape of the curve is very similar. This is no surprise, since pure neutron matter is 
like a Fermi gas, with attractive forces between neutrons lowering the energies with respect 
to the non-interacting system. The search for any kind of phase transition in the range of 
5 — 50 MeV was to no avail. It is likely that a phase transition occurs at lower temperature. 
We finally calculate two additional observables and compare them with other calculations 
found in the literature. The symmetry energy, which appears as a coefficient agym in the 
semi-empirical mass formula 

^^y^ _ {N-Zf 

is plotted in Fig. ^ We used the energy per particle of pure neutron matter (PNM) and 
SNM, subtracted them and interpolated the result on a mesh. Since the error bars for pure 
neutron matter are larger for low temperatures, the graphs should be viewed with caution. 
Nevertheless it appears that the symmetry energy is increasing with density and decreasing 
with temperature, as one would expect. Indeed, at high temperature SNM and PNM both 
are more like a Fermi gas, and only at low temperatures do they become different. The 
observed dependence on density can be explained by the fact that a dilute system is barely 
interacting while the probability of clustering increases with density. At saturation density 
and low temperature, we obtain a coefficient of 

a,ym ^ 38 ± 3 MeV, (54) 



which is not too different from the generally accepted value [p^] of a 



"sym 



28.1 MeV. This 



discrepancy is in part due to the fact that the calculations for pure neutron matter have not 
converged completely. 

The first sound velocity has been calculated using the formalism of relativistic fluid 
dynamics: 



u/c 



\ 



dp 
de 



(55) 



where 



p X rriNC + 



E 



(56) 



and 



P = P 



dE/A 



dp 



(57) 



In general, the first sound we obtain is too low compared to Ref. ||6|j28[|, but the velocities 
are of the same order of magnitude. Several calculations P,^ show a violation of causality 
at a few multiples of the saturation density, and so do ours. Our results (Fig. |^) show 
the correct temperature dependence in the sense that it conforms with the compressibility: 
higher sound speed for intermediate temperatures (high incompressibility) and lower speeds 
for both low and high T. 



IV. DISCUSSION AND CONCLUSION 

The model considered in this paper is an exact treatment taking first steps towards a 
more realistic Hamiltonian, and it is an improvement compared to previous calculations. 
Nevertheless it can be extended to include more physics, and details in the algorithm can 
be improved. First of all, more computer power is necessary to reduce finite size effects. A 
lattice of 10 x 10 x 10 points would be desirable, and also the imaginary time dimension could 
be pushed further. This requires stable matrix techniques. The present code can handle 30 
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time slices comfortably using commonly known sparse matrix techniques. But, as the lattice 
spacing decreases, one needs to go to larger imaginary time to separate the ground state 
from excited states. At the same time it is not possible to increase A/5 as it would induce 
finite time effects. Therefore, an improved effective matrix algorithm would be needed to 
allow for more matrices to be multiplied. Along with a bigger lattice, the resolution of the 
potential could be increased, resulting in next-to-nearest-neighbor and further interactions. 
This extension of the spatial dependence of the potential can easily be accomplished and is 
only restricted by computational power. 

Another big hurdle is the sign problem. The solution to this obstacle will result in a 



huge advancement in many areas of computational physics and chemistry. Rom et al. ||2^ 
have made some progress which could prove beneficial for the model described here too: It 
basically consists of shifting the contour of the auxiliary field integrals, which is equivalent 
to subtracting a mean-field from the problem. We plan to investigate this method and its 
application to nuclear matter more rigorously in the future. 

The physics of nuclear matter itself is certainly more involved than the current model can 
account for. Mesons are not included as explicit degree of freedom, and the various exchanges 
are only simulated indirectly through the choice of potential and its parameters, very much 
like in AV18 or other potentials. Realizing that the auxiliary fields behave like massless 
bosons, one could ponder how a Monte Carlo procedure would look like that includes meson 
exchange directly. Such a procedure could be quite similar to already established auxiliary 
field Monte Carlo procedures. 

It is known that three-body and perhaps higher-order many-body forces are important 
to describe saturation properties of nuclear matter correctly. However, incorporating these 
forces in a Monte Carlo calculation is currently impossible, basically because there is no 
scheme to reduce higher-order forces to the single particle formalism. Such a scheme could 
lie in a multiple application of the Hubbard-Stratonovitch transformation for a single many- 
body interaction. But as long as such a scheme is not available, an approximation could 
be established on top of this two-body calculation that incorporates higher-order effects. A 
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first attempt would be to calculate the three-body contribution to the energy obtained from 
the one-body densities of this Monte Carlo calculation and a given three-body Hamiltonian. 
In conclusion, this project has produced promising results which should be viewed as a 
starting point to an exact solution of infinite nuclear matter. In a model with a relatively 
simple Hamiltonian, and further limited by a very small lattice, we were able to reproduce 
saturation properties of symmetric nuclear matter. The energy of pure neutron matter, 
using the same potential, gave reasonable results, even though it had not yet converged. 
Furthermore, we presented evidence in form of mechanical and thermodynamical observables 
which support the existence of a phase transition from a Fermi gas to a clustered system. 
Particle fluctuations of the gas phase seem to reach a maximum at T ^ 14 MeV. The heat 
capacity and compressibility also have a maximum at around this temperature. Entropy 
and grand potential show a behaviour as it is expected for a first-order phase transition. 
Other quantities like symmetry energy and first sound velocity show reasonable agreement 
with other calculations. 
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FIGURES 
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FIG. 1. E/A for symmetric nuclear matter as a function of density p and for different temper- 
atures. The purpose of the hues is to guide the eye. 
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FIG. 2. Grand canonical potential of symmetric nuclear matter for different chemical potentials 
fi. The solid lines represent the potential for a noninteracting Fermi gas in continuous space. 
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FIG. 3. Density fluctuations in symmetric nuclear matter. The upper panel displays the mod- 
ified density / while the lower panel shows the derivatives of / with respect to chemical potential 
^ which are proportional to the fluctuations. 
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FIG. 4. Compressibility of symmetric nuclear matter. The minima of the energy curves have 
been fit to a parabola with a x^ ^ 1-5 per degree of freedom. 
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FIG. 5. Heat capacity and entropy for a finite piece of symmetrical nuclear matter. The two 
graphs on the left show the case n = 0.0 MeV, the right ones for fi = 4.0 MeV. The heat capacity 
(upper panels) shows a distinct maximum, the entropy (lower panels) a relatively sharp drop with 
decreasing temperatures. 
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FIG. 7. Symmetry energy for symmetric nuclear matter as a function of density and temper- 
ature. Shown is the coefficient asym of the semi-empirical mass formula. The left panel shows a 
contour plot, the right one shows one-dimensional cross-sections of it at different temperatures. 
asym is increasing with density and decreasing with temperature. 
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FIG. 8. First sound velocity for symmetric nuclear matter. The temperature dependence of 
the speed corresponds to the compressibilities as shown in Fig. Q. 
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